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A matrix method of analysis is developed for mildly nonlinear, 
multiple -input, multiple-output systems with memory (e.g., nonlinear 
multiport networks and multichannel communication systems). The 
method is based on a Volterra- series representation whose kernels 
are two-dimensional matrices rather than multidimensional arrays. 
This is made possible through the use of the Kronecker product of 
matrices, which results in a compact formulation. The response of 
the aforementioned systems to multiple sinusoidal excitations is also 
studied. Moreover, formulas are given for various system operations 
(e.g., addition, cascading, inversion, and feedback), which can be 
used to describe a complex system as an interconnection of simple 
subsystems. 

I. INTRODUCTION 

Communication, control, and instrumentation systems employ com- 
ponents, such as amplifiers and mixers, which are inherently nonlinear. 
Even when the nonlinearities are nuld, as is often the case, they can 
produce bothersome signal distortion that limits the system perform- 
ance. The nonlinear components themselves, and the other linear 
components used in the system, are generally frequency- dependent, 
i.e., they have memory. Numerous studies are available in the literature 
for the analysis of mildly nonlinear systems with memory through the 
use of Volterra-series expansions. 1 " 21 The classic paper by Bedrosian 
and Rice, 7 the recent paper by Chua and Ng, 14 and the book by Weiner 
and Spina 17 cover that subject very thoroughly. Also, the paper by 
Gopal, Njakhla, Singhal, and Vlach 12 is interesting in that it evaluates 
the range of accuracy of the Volterra-series approach by comparing it 
with a nearly exact, but quite involved, method of analysis. The book 18 
and paper 19 by Schetzen deal mainly with random inputs. The condi- 
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tions for the existence of a Volterra-series representation have recently 
been studied rigorously by Sandberg. 21 

For the most part, the studies mentioned above are limited to 
systems with one input and one output, i.e., "scalar" systems. This 
scalar representation is usually not easily applicable to Multiple-Input, 
Multiple-Output (mimo) systems. Such systems include, for example, 
nonlinear multiport networks, multichannel communication systems, 
and transmitting or receiving systems employing multibeam antennas. 
In principle, one can represent these systems by a set of dependent 
scalar Volterra equations. This was done, for example, in the papers 
by Narayanan 3 and by Bussgang, Ehrman, and Graham, 9 where node 
equations were used to analyze nonlinear, two-port network models of 
bipolar transistor amplifiers. This method of analysis is tractable only 
when the numbers of nodes and of nonlinear elements in the network 
are small. For example, when the above authors considered the analysis 
of two-stage transistor amplifiers, they were forced by the complexity 
of the cascade equations involved to assume that the interaction 
between the stages, i.e., the loading effect of one stage on the other, is 
linear. While this might have been a reasonable approximation in their 
particular case, it is not valid in general. A symbolic matrix inversion 
algorithm that simplifies the computational aspects of the nodal 
method of analysis was recently discussed by Thapar and Leon. 15,16 

To conveniently handle the problem of two-port networks, or to 
analyze nonlinear multiport networks in general, one needs to use a 
black-box representation of the network, as is usually done in linear 
networks. For example, consider a nonlinear, two-port network, which 
has two independent port variables (e.g., the port currents) and two 
dependent port variables (e.g., the port voltages). One should be able 
to express the latter variables in terms of the former (e.g., by a 
nonlinear impedance representation). Furthermore, one should be able 
to perform transformations among various network representations 
(e.g., from impedance to cascade parameters), and to carry out the 
computations involved in interconnecting several networks together to 
form a complex network (e.g., through cascading). The same operations 
are also needed in the analysis of other nonlinear mimo systems. 

The purpose of this paper is to develop a method for analyzing 
mildly nonlinear mimo systems with memory. This method, which 
employs Volterra-series whose kernels are two-dimensional matrices, 
facilitates the systematic performance of various useful system oper- 
ations, such as addition, cascading, inversion, and feedback. The 
application of the results of this study to the analysis of mildly 
nonlinear multiport networks will be the subject of a future paper. 

Actually, Weiner and Naditch, 10 and Gopal, Nakhla, Singhal, and 
Vlach 12 used multidimensional arrays of Volterra kernels to represent 
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nonlinear, two-port networks. The same was suggested by Chua and 
Ng 14 for extending their results to multiple-input systems. All of these 
analyses can also be generalized to multiport networks and other mimo 
systems. The resulting notation is similar to the index notation dis- 
cussed in the beginning of the next section and in Appendix A. This 
notation, though more natural in its initial formulation, turns out to 
be cumbersome when attempting to perform the aforementioned sys- 
tem operations. 

II. REPRESENTATION OF NONLINEAR MEMORYLESS MIMO SYSTEMS 

A nonlinear, memoryless scaler system is characterized by its in- 
stantaneous input-output transfer function. When this function is 
analytic, as is the usual case encountered in practice, it can be repre- 
sented by the power-series expansion 



w = P (U « + P m u 2 + P ( V + 



(1) 



where u = u{t) is the input, w = w{t) is the output, and P (k) , k = 1, 2, 
3, • ■ • , are system constants. The corresponding representation of a 
nonlinear, memoryless, mimo system with n inputs, Uj = uj(t), j = 1, 2, 
■ ■ ■ , n, and m outputs, w, = Wi(t), i = 1, 2, ■ • • , m, is 



+ 1 1 1 Pffih&hH&k + ■ ■ • , i = 1, 2, 
/i-i /2=i kr 1 



, m, 



(2) 



where P*!..j», k = 1, 2, 3, ■ ■ ■ , are (k + l)-dimensional m x n x 
x n arrays of system constants. The notation used in (2) will be 
referred to as the "index notation." It is similar to that used in Refs. 
10 and 12, but the superscripts and subscripts are interchanged. We 
now proceed to represent (2) in the "matrix notation." 
Let 



u = u(t) = 



Ui(t) 

u 2 (t) 



, w = w(f ) = 



Wi{t) 

w 2 (t) 



(3) 



Unit)] \w m {t)_ 

be the n X 1 and m x 1 input and output vectors, respectively. The 
first (i.e., linear) term in (2) can be written as an ordinary product of 
matrices in the form w = P (U -u, where P U) is the m X n matrix 
[Pi/*]- We will now show that the remaining terms in (2) can also be 
written in a matrix form through the use of the Kronecker product of 
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matrices. 22 " 24 Appendix B defines this product and gives some of its 
useful properties. Actually, Harper and Rugh 11 employed the Kro- 
necker product in conjunction with state variables to study factorable, 
scalar, nonlinear systems. Also, Brockett 25,26 used a reduced form of 
the Kronecker product (to be explained shortly) in the state-variable 
representation of scalar, time-varying, nonlinear systems that are 
linear in the control variable. 

As is explained below, the elements of the (k + 1) -dimensional, 
m X n X •-- X n arrays, {P£\..>}, can be reorganized to form 
two-dimensional, m X n k matrices, {P (A) }, such that (2) can be written 
in the matrix form 



w = P (,, 'U + P (2| '(uxu) +P (31 -(uxux u) + 



(4) 



where "X" is the Kronecker-product sign. As mentioned in Appendix 
B, we will employ left Kronecker products. 

To understand (4), we note from (61) that the A-fold Kronecker 
product u X u X • ■ • x u results in an n* X 1 vector whose j th element 
is given by 



where ji,^, 



[u x u X • • • X u] ; = u jl u j2 • - ■ u jh , 
, jk are uniquely determined from 



j=ji + n(j 2 - 1) + ••■ + n* Hjk - 1) 



(5) 



(6) 



Thus, to make (4) equivalent to (2), the i-j element of the m X n k 
matrix P (A) should be given by 



[P ik) h = P 



I/l •••/*! 



(7) 



where j is given by (6). 
For example, if m - n - 2, (5)-(7) give 











U1U1U1 








«2«1«1 








~U\U~ 




«1«2«1 






uxu = 


«2«1 
«1«2 
«2«2 


,U X U X u = 


U2l'.2U\ 
«1«1«2 
U2U1U2 


> 


(8) 




U1U2U2 








«2«2«2_ 






p«2) _ 




5 


(9) 


p(3> _ 


P&v 


PWnP?&P% 
P^iP^iPl 


!2lPlll2 
!21~2112 


Pi^Pi 

p^pi 


3) p(3> 
122T1222 
3) p(3) 
122*2222 


(10) 
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Note that each of the Kronecker- pro duct vectors given in (8) has 
redundant entries. In Brockett's notation cited earlier, these entries 
would be removed. For example, the 4x1 vector, u X u, would be 
replaced by the 3x1 vector [«?, U1U2, ut], and the corresponding 2 X 
4 matrix, P (2) , given in (9) would be reduced to a 2 x 3 matrix, etc. 
However, when the system has memory, no redundant entries occur, 
since, as can be seen for example from (11), one needs to evaluate 
Kronecker products of the form u(t x ) X u{t 2 ), etc., where *i ¥> t 2 . 

In the remainder of the paper, we will employ the compact matrix 
notation used in (4) rather than the index notation used in (2). 
However, on some occasions, it is helpful to keep track of the interre- 
lation between the two notations. Thus, some key equations in the 
paper are rewritten in Appendix A in the index notation. 

III. REPRESENTATION OF NONLINEAR MIMO SYSTEMS WITH MEMORY 

The usual Volterra-series expansion used to represent nonlinear, 
time-invariant, scalar systems with memory 1 " 21 can be generalized 
through the use of the notation of (4) to represent mimo systems by 
the matrix equation 

w(t)= p (1) (ti).u(£-ti)c?Ti 

+ P <2) <T,, T 2 )-[U(^ - Tl )XH(f- T 2 )]rfT!rfT 2 



\\\. 



P (3) (T1, T2, T 3 ) 
[U(t — Ti) X U{t — T 2 ) X u(t — T 3 )]dTidT 2 dT3 



(ID 



where u(t) and w(t) are, respectively, the n x 1 and m X 1 input and 
output vectors given by (3), and where p (A) (Ti, • • • , r k ), k = 1, 2, 3, 
• ■ ■ , are two-dimensional, m X n k matrices of system kernels. Note 
that if P ( *>(ti, ... , ta) = P^Sin) ... d(r k ), where 8(t) is the unit 
impulse function, then the system becomes memoryless, and (11) 
reduces to (4). 

As is the case for linear systems, it is more convenient to represent 
(11) in the frequency domain. To do this, we introduce the dummy 
time variables, ti, fe, • ■ - , ?*, and rewrite the k th order output com- 
ponent in (11) as 
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•[u(*i ~ Ti) X - • ■ X U{t k - T A )]dTi ■ • ■ dr k . (12) + 

Thus, (11) becomes 

W (*) = w (1) (i) + w (2, <*, + w»ft t, *)+•■• . (13) 

Now, we introduce the single-dimensional Fourier-transform pair 

X(/)=| x(t)expi-j2vft)dt, (14a) 

*/ — OD 

x(t)-[ X(f)exp(j^ft)df, d4b) 

to represent the transformations u(0 <-* W) and w(#) «-> W(/). 
Similarly, we introduce the multi-dimensional Fourier-transform pair 

Y(/i, ■••,/*)-( ••■ J yte> ■••.*) 

•exp[->27r(/i^i + • ■ ■ + fktk)]dti ■ - - c#a, (15a) 
y(*i, ■■-,«*)=[ ■■• J Y(/i, -..,/*) 

■ exp[y2ir(M + ■ ■ ■ + /*fa)]d/i ■ ■ • djf*. d 5b ) 

to represent the transformations p (A) (ti, • • • , t*) <-> P (A) (/i, • • • , /*) 
and w {k \t u • • • , t h ) ** W (A) ( A, ■ ■ • , A). It can be shown from (14) and 
(15) that (12) can be written in the frequency domain as (see Refs. 1 
through 4, 7 through 9, 14 and 17) 

W»»</ lt .-.,/*) =P W ( A, ■•-,A)-[U(/i)x ••■ xU(A)]. (16)+ 
The Fourier transform of the output becomes 



W(fl = W {1) (f) + W< 2) ( A, f - AMA 

J — oo 



I 



J 



f W (3 >( A, A, f- A - A)^ArfA + • ■ ■ - (17) 



+ All equations in the paper marked by a dagger are rewritten in the index notation 
in Appendix A, where the same equation numbers are used. 
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Fig. 1 — A nonlinear, time-invariant, mimo system with memory having n inputs and 
m outputs. 



(ID 




(14} 



Fig, 2 — Interrelations among the input, output components, and total output in the 
time and frequency domains for the nonlinear mimo system of Fig. 1. The numbers in 
parentheses represent equation numbers in the text. 



Note from (13) and (17) that the single-dimensional Fourier transform 
of w lk) (t, ■ • • , t) is given by the &th term in (17), which is not equal to 
W (A, (/-, ... , /) unless k = 1. 

A schematic diagram of the system represented by (11) through (17) 
is given in Fig. 1. The interrelations among the input, the output 
components, and the total output in the time and frequency domains, 
and the corresponding equation numbers, are indicated in the flow- 
chart of Fig. 2. 

IV. KERNEL SYMMETRIZATION 

The representation of the response of a nonlinear scalar system to 
sinusoidal and Gaussian excitations is greatly simplified if each of the 
kernels, P {k) /i , • • • , f k ), or equivalently, p m (ti , • ■ ■ , t* ), is a symmetric 
function of its arguments. 7 " 9,14 The generalization of this symmetry 
requirement to nonlinear mimo systems is somewhat more involved. 
Following the reasoning given in the aforementioned references, one 
can show that it is the output components, W (A) (/i, > • ■ , f k ) given by 
(16), or equivalently, w'*^^, • ■ • , tt) given by (12), that are required 
to be symmetric functions of their arguments. For example, for k = 2, 
it is required that W (2) (/i, f 2 ) = W (2) (/ 2 , f\)\ and thus, from (16), 



P'*(/i, £)■(!!, X Us) = F 2, (/2, /l)-(U2 X u,), 



(18a) 



where U(/i) is replaced by u, for generality. Similarly, for k — 3, it is 
required that W (3) (/i, f it f 3 ) = W m (f a , fa f r )l and thus, from (16), 
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T? m (fu fz, f*)-(*i x U2 x u 3 ) = P™(f a , k, f y )-(u a XUpX u y ), (18b) 

where a, ft, y assume all permutations of 1, 2, 3. For a scalar system, 
(18a) and {18b) are indeed equivalent to requiring the corresponding 
system kernels to be symmetric functions of their arguments, as 
mentioned above. 

To find the symmetry requirement implied by (18) on the kernels of 
a mimo system, we need to introduce the n 2 X n 2 "reversing" matrix, 
R, and the six n 3 x n 3 "permutation" matrices, O a ^ y , where a, ft y 
assume all permutations of 1, 2, 3. These matrices have properties such 
that if ui , u 2 and u 3 are n X 1 vectors, then 

R-(uiXu 2 ) = 112X11!, (19) 

® a (iy - (ui x u 2 x u 3 ) - u„ x ii/, x u T . (20) 

Appendix C defines these matrices and gives some of their useful 
properties. 

Finally, (18) through (20) give the required symmetry conditions of 
the kernels as 

P (2, (A,/2)=P (2, (/2,/i)-R, (21a) + 

P i3} (fu k, *> -P»</« U, f?)'*afly. (21°)+ 

The generalization of (21) to higher-order kernels requires the intro- 
duction of permutation matrices of more than three jndices. 

If the given system kernels, say, P (2) (A, ^) and P (3) (A, £, fs), are 
unsymmetric, they can be symmetrized, i.e., made to satisfy (21), 
through the use of the relations 

P (2, (/i, / 2 ) = fc \P*(fu h) + P l2) (f 2 , /5).R], (22a) + 

P (3, (A,A,A) = , / 6 2 P™(f at £,£>•*«*, (22b) + 

where the summation is performed over a, /?, y assuming all 6 permu- 
tations of 1, 2, 3. These symmetrization relations are generalizations 
of those discussed in Refs. 7, 9, and 14 for scalar kernels. 

V. RESPONSE TO SINUSOIDAL EXCITATION 

The response of a nonlinear scalar system to multiple-sinusoidal 
excitation has been studied by several authors including Bedrosian 
and Rice, 7 Goldman, 8 and Chua and Ng. 14 Here we generalize some of 
their results to nonlinear mimo systems. 

5.1 Multiple-exponential excitation 

Let the input vector be 
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u<0= S u ( -exp(y27r/rf), (23) 

i-i 

where the u/s are time-independent, complex, n X 1 vectors. The 
Fourier transform of u(t) is 

U(/) = Su;S(/-£). (24) 

i-l 

Substituting (24) into (16), and using (15b), one obtains the kth order 
output component 

w m (t) = w^u, • ■ ■ , t) = I ■ ■ - 2 (CP W (A. ■ ■ ■ , /y 

'i-l «*-! 

■ (Ui, X ... X u ik )]exp[j2<jT(fi l + • • ■ + 4)*]}. (25) 

Finally, the output, w(t), is obtained from (13), i.e., by summing w {h) (t) 
from k = 1 up to any desired order. Note that (25) is valid whether 
or not the system kernels are symmetric. 

5.2 Single-frequency excitation 

Let the (real) n X 1 input vector be 

u(t) = Real[a exp(j2irft)] 

- & a exp(j2irft) + Msa*exp(-/27r#), (26) 

where the asterisk refers to complex conjugation. Comparing (26) to 
(23), one obtains / = 2, Ui = x h a, u 2 = V2 a*, f\ = f, and ft = —f. Thus, 
using (25), and assuming that the kernels are symmetric, i.e., that (18) 
is satisfied, one obtains the following expressions for the various Ath 
order output components: 

w' 1 ^) = % [P iu (f).a]exp(j27rft) 

+ K [P {1) (-f).a*]exp(-j2irft). (27a) 

w {2) (t) = x h [P m (f, -/)-(a x a*)] *- (d- c term) 

+ V 4 [P (2) (/, /)-(a X a)]expU27r(2f)t] 

+ %[P m (-f, -/)-<a* X a*)]exp[-y27r(2/)^]. (27b) 

w W(t) = % [P <3) (/", f, -f)-(a. x a X a*)]exp(j27rft) 

+ % [P (3) (-f, -f, /)-(a* x a* x a)]exp(-/27r#) 

+ V« [P (3) (/; /, 0-(a x a x a)]exp[j27r(3nt] 

+ Va [P (3) (-/, -/, -f)-(a* X a* x a*)]exp[->27r(3A)^. (27c) 

Note that the asterisks on the a's correspond in number and location 
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to the negative signs in the frequency arguments of the associated 
kernels. 

If the system is real, i.e., if P (A) (ti, ■ • • , t*), k - 1, 2, 3, • • * , are real, 
then it can be shown from (15) that 

P w (/i, •■■,/*> = [P <A) (-/i, - • * , -A)]*- (28) 

As expected, (28) implies that all the output components given in (27) 
are real. In that case, it can be shown through generalizing (27) that 
the total mth harmonic output term, w m (£), m - 0, 1, 2, ■ ■ • , is given 

by 

k 



w m (() = e m Real 



exp(y'27rm/*) £ 



,_ A i k - m 



(29) 



■P (A) ( f • ■ • , f,-f,"->-f)M {k+m)m X <a*) [<A - m,/2J ] 

ik+m)/2 (k-m)/2 ) 

where a l/] is the /-fold Kronecker product ax ... X a, € m is the 
Neuman factor (which is equal to 1 when m = 0, and is equal to 2 
when m/0), and P (0) is defined to be zero. 

Because of the symmetry conditions of (21), the kernels used in (27) 
and (29) satisfy the relations 

P ,z) (/ 1 /) = P ,2) (f,/)-R- (30a) + 

P <2) (/;-/ r ) = [P <2) (/,-/)]*-R. (30b) + 

p <3) (/; /, -/) = p (3, ( f, f, -n-*2i 3 , (30c) + 

P (3) ( f, f, f) = P (3) ( /, f, f) •*•»■ (30d) + 

In addition to the kernel symmetry requirement, (30b) is based on the 
assumption that the system is real, i.e., that (28) is satisfied. The 
implication of (30) is that the elements of each of the system kernels 
are not all independent. For example, if n = 2, (30a) through (30d) 
imply, respectively, that (i) columns 2 and 3 of P (2) ( f, f) are equal; (ii) 
column 2 of P (2) (/, —f) is the complex conjugate of column 3, and 
columns 1 and 4 are real; (Hi) columns 2 and 3 of P <3, ( /", /, -f) are 
equal, and so are columns 6 and 7; and (iv) columns 2, 3, and 5 of 
P t3> ( ft f> f) axe equal, and so are columns 4, 6, and 7. It is worth 
mentioning that (30a) and (30d), respectively, would also be satisfied 
by P (2> and P <3) of the memoryless system represented by (4). 

5.3 Two-frequency excitation 

Let the (real) nXl input vector be 

u(t) = Real[a exp(j^U) + b exp(j27rf b t)]. (31) 
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We assume that the system is real, and that the kernels are symmetric, 
i.e., that (28) and (18) are satisfied. One can use (25) to obtain the 
output corresponding to (32) by following the same steps used to derive 
(29). The leading terms at some of the various output frequencies are: 

w(t) \ d - c a fc[P ,2, < f a , -fa) • (a x a*) + P< 2 >( f bl -f b ) . (b X b*)]. (32a) 
w(Ok ~ Real{expO'277/ a *)[P n) (/a)-a 

+ 3 /4P (3) ( fa, f a , -U) ■ (a X a x a*) 
+ %P (3 >(/ a , f b , -f b ) -(axbx b*)]}. (32b) 

wU)|2/ -/ 4 a Real{ 3 /4 exp[j2w(2/ a - f b )t] 

.P (3| (/J„,-/ t ).(axaxb')}. (32c) 

w(t)\ lfn±mfh * Real J 2- (/+m - n ( l + m J exp[j2ir(lf a ± mf b )i] 

■P' ,+ ""(/», ■•-,/ n , ±A, ■ : - ,±/ 6 ) -[a 1/1 x (b*) 1 " 1 ] J, (32d) 
where £ > and m > 0, and where we defined b + sb and b~ = b*. 



5.4 Three-frequency excitation 

Let the (real) nxl input vector be 

u{t) = Realfa exp(j2irf„t) + b exp(j2irf b t) + c ex-p(j2irf c t)]. (33) 

Again, we assume that the system is real, and that the kernels are 
symmetric. Following the same steps used to derive (29) and (32), one 
can obtain the following leading terms at some of the various output 
frequencies: 

w(0U-c « M*P ,2 U, -/.)-(a X a*) + P m (f*>, -ft) 

• (b x b*) + P m (f c , -f e ).(c X c*)]. (34a) 

wtf)k « Real{exp(>2irAO[P ,n (A)-a + %P m (fa, f a , -f a ) 

.(a x a x a*) + %P ,3) (/L- A, -ft) -(a xbxb*) 

+ %P (3) ( £, £, -/J ■ (a X c x c*)]}. (34b) 

wU) |r 8+ /i-f e a Real{% exp[>27r( f a + f b - f c )t] 

■P (3) (/«, ft, -fc)-(a X b x c*)}. (34c) 
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W(t)\kfj±lf b ®mf c 

* Real ( 2 -'* +f+ '"- 1 » [k t,l + , m)! ewUMkfa ± lf b ®mf c )t] 
k\l\m\ w 

k I m 

.[a w x (b*) [/1 x (c©)' m] ]j, (34d) 

where k,l,m> 0, but at least one of them being nonzero, and where 
the sign symbols ± and (±) are each consistent throughout the equation, 
but are otherwise independent. 

VI. SYSTEM OPERATIONS 

6. 1 Operational notation 

Let the input- output relations given in (11) through (17) be written 
symbolically as 

W m = {P£UoU„, (35) 

where "o" means "operating on." The frquency dependence has been 
omitted for simplicity. The subscripts n and m are included to empha- 
size the numbers of inputs and outputs. On some occasions, these 
subscripts will be eliminated. 

If the system is linear, i.e., if P m = for k > 1, the operation in (35) 
reduces to an ordinary matrix product. Thus, 

W= {P (1, }oU = P n) .U. (36) 

The operational notation of (35), and the three system operations of 
addition, cascading, and inversion, which are discussed in the next 
three subsections, form an algebraic structure that permits a shorthand 
description of complex interconnections of nonlinear mimo systems. 
The laws of this algebra 1 are identical to the algebra of linear systems 
(i.e., the algebra of matrices) with two important exceptions — the left 
distributive law does not hold, and the laws of multiplication by a 
scalar constant are more complex. 

6.2 Addition 

Two systems, (P^j and {Q£U, having the same number of inputs, 
n, and the same number of outputs, m, are said to be "added" if they 
share the same input vector, U„, and if their respective outputs are 
added to form the final output vector, W m . This operation, which is 
shown schematically in Fig. 3, is represented by 

2232 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1 982 



pi*) 



s> {<] 



=c> 



I V {C] 



Fig. 3— A schematic representation of the addition operation {&£„} = {PSi.} + 
{QtfiJ. 



W m = {P&}oU„ + {Q^„}oU„ 

= [{Fffi.} + {Qtfix}]oU„ 
= {S£UoU„. 
The kernels of the sum system, 

{SffiJ = {P<& } + {Q&}, 
are given by 

8 1 * 1 !/,, ...,/ A )-P w (^ •.-,/*) +Q w (/i. 



A), 



(37) 



(38) 



(39) 



where the plus sign refers to matrix addition. 

One can define a subtraction operation in an obvious manner. A 
multiplication operation, 1,14 which is more involved, can also be de- 
fined. 

6.3 Cascading 

When the output vector, W«, of a system, {PJ^,} , is used as an input 
vector to a second system, {Qi,m}, whose output vector is X/, the two 
systems are said to be in "cascade." This operation, which is shown 
schematically in Fig. 4, is represented by 



r 



^ p, 



' m.n 



,(*) 



T <*) 
V, n 



i * 



Fig. 4— A schematic representation of the cascade operation {!$} = {Q&} * (P^i.}- 

MATRIX ANALYSIS 2233 



= {QfS}o[{PS3.}oUj 

= [{Qffi} * {P^)]oU„ 

= {T}*>}oU„, (40) 

where the asterisk refers to the cascade operation. The kernels of the 
cascade system, 

{Tft»} = {Qffi} * {Pffi.}, (41) 

can be obtained by substituting the output expression of the first 
system into the system equations of the second system, as was done in 
Refs. 3 and 9 to derive the cascade relations of scalar systems. This 
procedure is straightforward, but somewhat tedious. A simpler ap- 
proach is to employ the harmonic probing method discussed in Refs. 
7 and 9, and the expression for the response of nonlinear vector systems 
to multiple-exponential excitation given in (25). The resulting relations 
for the cascade kernels are 

T (1) (A) = Q a, (A)-P (1) (A), (42a) + 

T< 2) (A, A) = Q (1, <A + A)-P <2) (A,A) 

+ Q (2) ( fa A) ■ [P m ( A) x P m (f z )l (42b)* 
T (3) (A, A, ft) = Q (1) (A + A + A>-P (3) (A, A, A) 

+ Q <2) (A, h + A)-[P (1) (A) x P< 2) (A, ft)] 
+ Q (2) ( A + fa A)-[P (2) ( A, A) x P (1) ( A)] 
+ Q <3) (A, fa A)-[P U) (A) x P (1, (A) x P (1 >(A)]- (42c) + 

A generalization of (42) for arbitrary k is given in Appendix D. 

If the kernels of the cascaded systems are symmetric, i.e., satisfy 
(21), then it can be shown that the resulting second-order kernel given 
by (42b) is also symmetric. However, the resulting third-order kernel 
given by (42c) is not symmetric. This fact is indicated by the presence 
of the circumflexes. 

As mentioned in Section IV, it is desirable to deal with symmetric 
kernels. Thus, using the symmetrization relation given in (22b), assum- 
ing that the kernels of the cascaded systems are symmetric, and 
employing the properties of the reversing and permutation matrices 
given in Appendix C, one obtains the symmetric form of (42c) as 

T< 3, (A, fa A) =Q (1, (A + h + A)-P <3) (A, fa h) 

+ 2 /3{Q (2) ( fa A + A)-[P a) ( A) x P°»( fa A)] 

+ Q <2) ( fa A + A)-[P (1) ( A) x P (2) ( A, A)]-*23i 

+ Q (2) ( A + fa A) ■ [P <2> ( A, A) x P a) ( A)]} 

+ Q (3) (A, fa A)-[P U) ( A) x P (I) ( A) x P (,) ( A)], (42c) t 

where <&23i is defined in (68) and (69). 
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If the first system, {Pm, n }, is Linear, (42) reduces to 



T lk) (fu .■■,£> = Q (A) (/i, ..-,A)-[P 0) (/i)x...xP (1) (/i)]. (43) 

On the other hand, if the second system, {Qi^}, is linear, (42) reduces 
to 



T (A >( /!,..- , A) =Q (A, (/i + •■• +/*).P (A) (/i, ••• ,/i) 



(44) 



6.4 Inversion 

Let the numbers of inputs and outputs in the system represented by 
(35) be equal, i.e., m = n. Suppose that it is required to find the input 
vector, U„, in terms of the output vector, W«. This inversion operation 
is represented by 

U„ = {PISJ-W. = {Qffl}oW„. (45) 

To find the kernels of the inverse system, 

{Qffi} = {F8r\ (46) 

it is helpful to use the interpretation given in Fig. 5, which defines the 
inversion operation in terms of the cascade operation and the identity 
system, {1„}, where 1„ is the n X n identity matrix. Thus, applying the 
symmetric cascade relations of (42) to Fig. 5b by interchanging the 
roles of P and Q, setting T (1) = !„, and T'* 1 = for k > 1, and solving 
for Q <A) , one obtains the symmetric inversion relations 

Q (1, (/i) = [P u, </i>r\ (47a) 



I 



o {«£) -(«} 



■N 



(a) 



r 



» {«} - {«}-' 



~v 



. 1. 



[■a] 



L. 



1=** 
K 



(b) 



Fig. 5— Two equivalent interpretations of the inversion operation {Qi^} = {P«* n } ': 
(a) (OKI) * CPJSi) = {In}, and (b) {Pi**} * {Qjft} = {1„}, where fl»} is the identity 

system. 
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Q (2) (A, h) = -Q (1, (/i +/ r a )-P w (/ r i, / 2 )-[Q (1, (/i) x Q (1, (/£)], (47b) 
Q <3) (/i, fa h) - -Q U) (/i +& +f 3 ).(%{P {2) (fi, f2+f 3 ) 
•[Q a) (A)xQ< 2) </ 2 ,/ 3 )] 

+ P< 2 >(/ 2 , f a + /,).[Q (1) (/ 2 ) X Q (2 >( fa A)]-«28i 

+ P (2) (A + fa / 3 ).[Q <2) (/i, h) x Q (1) (/ 3 )]} 
+ P <3, <A, fa / 8 )-[Q (1) (/i) x Q (lt </2> x Q (1 >(f 3 )]). (47c) 
Note that the inverse system exists if and only if P (1) ( f) is nonsingular. 

6.5 Feedback 

As an application of the three system operations discussed in the 
previous subsections, consider the nonlinear, feedback, mimo system 
shown schematically in Fig, 6, where both the forward, {PUJ,}, and 
reverse, {Qi*™}, branches are nonlinear. Using the operational notation 
of (35), one obtains 



X„ = U„+{Q£UoW m) 



(48a) 
(48b) 



where U„, X„ and W m are the n x 1 input vector, the n x 1 intermediate 
vector, and the m X 1 output vector, respectively. Substituting W m 
from (48a) into (48b), solving for X„ in terms of U„, and substituting 
the result in (48a), one obtains the feedback system equation 



W m = {FXMoUn, 



where 



|(A> 



{F«J.> = {P£U * [{!„} - {Q&} * {P™ }] 



><*) n-i 



(49) 



(50) 



Thus, the kernels of the feedback system can be obtained by applying 
the subtraction, cascade, and inversion operations discussed above. 
However, the explicit formulas for these kernels wul not be given here. 



{■»} 



\ n/n J 



l_. 



\ r m,n J 



Fig. 6 — A schematic representation of a nonlinear mimo feedback system. 
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Actually, special cases of these formulas have been obtained for scalar 
systems in Refs. 2, 6, and 7. 

Note from (50) and Section 6.4, that the feedback kernels exist if 
and only if the n x n matrix [1„ - QjJJU/") ■ Pm,U/)] is nonsingular. 
Note also that if m = n, and if P$,( /) is nonsingular, then (50) reduces 
to 

{Fsa-HPSr^cos&r 1 . ( 51 > 

If the system in Fig. 6 is changed to a negative feedback system, 
then the minus signs in (50) and (51) should be changed to plus signs. 

Vlt. CONCLUSIONS 

A method of analysis has been presented for mildly nonlinear mimo 
systems with memory. The method utilizes Volterra series whose 
kernels are two-dimensional matrices. The analysis was made possible 
through the use of the Kronecker product of matrices, which is a 
simple but powerful tool in matrix theory. This results in a compact 
representation of the system equations, and facilitates the systematic 
performance of various useful system operations, such as addition, 
cascading, inversion, and feedback. These operations can be used to 
describe a complex, nonlinear mimo system as an interconnection of 
simple subsystems. 
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APPENDIX A 
Index Notation 

Here we rewrite, in the index notation, some of the key equations 
marked by a dagger ( + ) in the body of the paper. The same equation 
numbers are used here as are used in the text. Before doing so, 
however, we note from (7) that, for mimo systems with memory, the 
matrix kernels used in the matrix notations are related to the array 
kernels used in the index notation by the relations 

[P w (ti, ■ - ■ , r*)]y = pf v . . Jh (r u ■ . ■ , t*), (52) 

[P (At ( h - • • , fk)lj = /f. . *< A, - - ■ , fk), (53) 

where y is given by (6). 
A list of the equations in question follows. 
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y"i-i hr 1 U-« J 

. u>,(*i - ti) • ■ ■ u h (t k - r k )dri ■ ■ ■ dTk \ (12) f 



•EM /l> ■■• LT A (/i). (16) + 

Pj?A*< A' A. A> = P ^W A, ft, fr>- ( 21b > + 

Pjf*< A /a = ^"a< a. a> + *S*< A. /»]■ (22a)+ 

P5**( A. ^ A) - 5 I ^w( A. fc A). < 22b > + 

p™ J2 (f,n = pZsf>n- < 30a > + 

P*MJLf.f.-f)=P?lj*if.f.-f). (30c) + 

rg(A)- 2 Qff(A)PS(/i). <42a) + 

a = ] 
m 

rSfc< A, = S Q?< A + /»F5*( A, A) 



f&tt</i. A. A) = S <?!"< A + A + f^u,Afu A, A) 

m m 

+ 12 [Q5<A.A + )»p$<A)p8a<A.A> 
+ QS(A + A, flPSUUi. A)P$<A)] 

m m m 

+ 2 £ S Q2i r (A,/ r 2,/i)Pi 1 J 1 (/i)PS(A)i > S(A)- (42c ) + 

tr-l 0-1 T =l 
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Tijjjz/sifh h\ h) — £ Qiaifl + ft + fa)Pmhh(fii h> fa) 

0=1 

+ .H [QS(A. * + h)P™mP [ $*jAh> h) 

+ QS(/2, h + A>p3(£)P£a<A, h) 

+ Qfttf, /i + f2)P^(fs)P^ j2 (fu m 






+ 11 I QZrifi* f*> fs)P t &(fi)P&(foP§(f*l (42c) f 

a-l £=1 y=l 



APPENDIX B 

Kronecker Product of Matrices 

Here we define the Kronecker product of matrices and summarize 

some of its properties that are used in this paper. More extensive 
coverage of this topic is given in Refs. 22 through 24. 

Let A = [ot > ] and B = [fc^/J be m a X n a and rrib X rib matrices, 
respectively. Their Kronecker product results in the m a mb X n a rib 
matrix, C = [c, cjf J, given by 



C= AxB= 



Ab u A612 
A 621 A622 



A6„ 



Ao„ 



A6i„„ 
Ab 2ni , 



Ab r 



(54) 



where "X" is the Kronecker-product symbol. Thus, 

where 

i c = i a + m a (ib — 1) , 
jc - y a + n a (./6 - 1). 



(55a) 

(55b) 
(65c) 



Note that, since i a ^ m a and ,/ < ra a , (55b) and (55c) have unique 
solutions for i a , ia , / a and jb in terms of j c and j c . Actually, (54) and 
(55) define the left Kronecker product. 22 One can also define a right 
Kronecker product, 23 * 24 which, however, is not used in this paper. In 
general, 



AxB/BxA. 



(56) 



It can be shown that the Kronecker product has the following 
properties: 



2240 THE BELL SYSTEM TECHNICAL JOURNAL, NOVEMBER 1 982 



A x (B x C) = (A x B) x C = A x B x C. (57) 

(A + B) x C = (A x C) + <B x C). (58) 

A x (B + C) = (A x B) + (A x C). (59) 

(A-B) x (C.D) = (A x C)-(B x D). (60) 

{A x B)" 1 = A" 1 x B 1 . (61) 

(AxB) r = A T xB T (62) 

In the above equations, "T" refers to matrix transposition, and the dot 
implies ordinary matrix multiplication. The dimensions of the various 
matrices are arbitrary, but of course, should be consistent with the 
requirements of the inversion, addition, and ordinary multiplication 
operations, where applicable. 

APPENDIX C 

Reversing and Permutation Matrices 

Here we define the n 2 x n 2 reversing matrix, R (n) , and the six n 3 x 
n 3 permutation matrices, 0$ T , which satisfy (19) and (20). The super- 
script "(n)" is used in this appendix to emphasize the dimensions. It 
can be shown from (19) and (55) that R m is given by (cf. Ref. 24) 

B2Lu-i>.**ui-ii = M». i> 7, *, / = 1, 2, • • • , n, (63) 

where S u „ is the Kronecker delta, which is equal to 1 if a = /?, and if 
a^jfl. For example, 

1000 

0010 

100 

0001 



R <2) = 



(64) 



It can be verified that 

R<»> = [R<"i] r =[R<">]- 1 , (65) 

where "T" refers to matrix transposition. Moreover, if Mi and M 2 are 
m x n matrices, then 

R (m> . (Mi x Ma) - R M = M 2 x Mi , (66) 

which is a generalization of (19). 
It can be shown from (19), (20) and (60) that 

4& = r*&] r = [0&]" 1 = U x R (n> , (67b) 

•ft = [*^] r = [^] _1 = R (n) X 1», (67c) 
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*& = [*S?i] r = [*Si] _1 = [U x R< n) ].[R (n) x l n ], (67d) 

*!Si = [*^] T = [*&i] _1 = [R (n) x l„]-[l n x R (rt >], (67e) 

•& = [*iSl] r = t*^!]" 1 = [l n x R ( ">] . [R<"> x 1„] . [l„ x R iny ] 

= [R (n) x l„].[l n x R (n) ].[R (n) x 1„], (67f) 

where l n and l„a are the n X n and n 3 x n 3 identity matrices, 
respectively. Also, it can be shown from (20) and (55) that, if a, /?, y are 
any permutation of 1, 2, 3, then the i-j element of <&$ r is given by 

where 

i = i a + n{i $ - 1) + n 2 (i y - 1), (68b) 

/ = /i + "(72 - 1) + »*(/a - 1). (68c) 

and where U , j\ , i 2 , j% t is, j* « 1, 2, ■ ■ • , «. For example, (67d) and (68) 
gives 

1000000 
00100000 
00001000 
00000010 
01000000 
10 
00000100 
L0 lj 

It can be verified that if Mi , M 2 and M 3 are mX n matrices, then 

*^-(M! x M 2 x M3)-[*^ r ] r = M a xM^x ML,, (70) 

which is a generalization of (55). Also, if M and N are m X n and 
m 2 X n 2 matrices, respectively, then 

•$ ■ (M x N) - [*£\] T = N x M, (71a) 

*S$.(N X M).[*M] r = M x N. (71b) 

Moreover, if M and K are m x n and m X n 2 matrices, respectively, 
then 

R (m) .(M x K)-[*&] r = K x M, (72a) 

R (m > • (K x M) ■ [#M] r = M X K . (72b) 

Finally, if M and L are m X n and m 2 x n matrices, respectively, then 

tt&MM x L).R< n) = L x M, (73a) 

*^-(L x M)-R (n) = M x L. (73b) 



*Ii = [*&] r = 



(69) 
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APPENDIX D 

General Cascade Relation 

Here we give a generalization of the unsymmetric cascade relations 
given in (42a) , {42b), and (42c) for arbitrary k, (cf. Refs. 7 and 9 for 
scalar systems) 

f^ifu ...,/*) = i I *i q (/) (/i + ■•■ + &, fa 

\ {h l +k 2 +- ■ -+ki=k) 

+ ..• + fk l+ k 2 , ■■-, fk-k l+ i + ■■■ + fk)V? {h Hfi, • • • , M 

XP^(/ ftl+1| .... f k]+ k,) X ... xP ikl) (fk-k l+ u ■■■ , A)] | (74) 

Note that the second summation contains f . _ J terms, and that the 

frequency arguments always appear in the order /i, f 2 , • • ■ , fk- As is 
the case with {42c), the cascade relation of {74) does not preserve 
kernel symmetry for k > 3. The symmetric form of (74), which would 
generalize (42c), will not be given since it requires the use of permu- 
tation matrices of more than three indices, which have not been 
introduced yet. 
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